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ABSTRACT 

Simulations estimating the brightness temperature (5Tb) of the redshifted 21 -cm from the 
epoch of reionization (EoR) often assume that the spin temperature (T s ) is decoupled from 
the background CMB temperature and is much larger than it, i.e., T s ^> Tomb- Although a 
valid assumption towards the latter stages of the reionization process, it does not necessarily 
hold at the earlier epochs. Violation of this assumption will lead to fluctuations in 5T\> that 
are neither driven by density fluctuations nor by Hn regions. Therefore, it is vital to calculate 
the spin temperature self-consistently by treating the Lya and collisional coupling of T s to 
the kinetic temperature, T^. In this paper we develop an extension to the BEARS algorithm, 
originally developed to model reionization history, to include these coupling effects. Here we 
simulate the effect in ionization and heating for three models in which the reionization is 
driven by stars, miniqsos or a mixture of both. We also perform a number of statistical tests to 
quantify the imprint of the self-consistent inclusion of the spin temperature decoupling from 
the CMB. We find that the evolution of the spin temperature has an impact on the measured 
signal specially at redshifts higher than 10 and such evolution should be taken into account 
when one attempts to interpret the observational data. 

Key words: Physical data processes: radiative transfer - cosmology: theory - observation - 
diffuse radiation - radio lines: general. 



1 INTRODUCTION 

Physical processes that occur during reionization are numerous 
and complex. Nevertheless, ionization of neutral gas (hydrogen & 
helium) and heating of the inter-galactic medium (IGM) can be 
considered the two primary influences of radiating objects during 
reionization. 

Currently, the most promising "direct" probe of reionization 
is the redshifted 21 -cm radiation emanating from neutral hydrogen 
during the epoch of reionization (EoR), whi ch are to be measured 
using upcoming telescopes like LOFAlfp] MWfl PAPElQand 
21CM/4j The intensity of the observed 21 -cm radiation depends 
on the ratio between the number density of electrons in the hyper- 
fine states in the ground state of a neutral hydrogen atom. This ratio 
is normally expressed in terms of the so-called 21 -cm spin temper- 
ature, T s . At the onset of the formation of the first reionizing ob- 
jects the spin temperature is equal to the CMB temperature since 
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at these redshifts the ratio between excited and ground hyperfine 
state electrons is completely determined by the CMB. However, as 
the number of ionizing sources increases, T s starts departing from 
Tcmb; slowly at the beginning, then rapidly approaching values 
larger than Tcmb . This evolution is typically ignored in most pre- 
vious studies of reionization which assumes T s ^> Tcmb at all 
times JMellema et al.|2006||Abel & Norman|2000| |2002l |Bromm| 
|et al.|2002l|Yoshida et al.|2003] >. 

Recently, Baek et al. (2009} have relaxed this assumption on 
T s at the dawn of reionization and explored its impact on the bright- 
ness temperature. They found a considerable considerable devia- 
tion from assuming T s ^> Tcmb at the beginning of reioniza- 
tion. Towards the end of reionization though, this assumption holds 
ground. But, in order to track the evolution of T s accurately, like in 
Baek et al. ( 2009}, it is necessary to perform a detailed 3-D Lya ra- 
diative transfer calculation. The Lya photons undergo a large num- 
ber (~ 10 5 ) of scatterings even in a marginally neutral medium 
before it is sufficiently off line-centre to "free stream". The scatter- 
ing angle after each encounter is completely random and therefore 
the radiative transfer is often done in a Monte Carlo sense < |Tasit-| 
siomi 2006, Baek et al. 2009, Pierleoni et al. 2009 ; |Laursen et al.| 
2009) to capture this random nature of Lya scatterings. 

Unfortunately these Monte Carlo radiative transfer schemes 
are computationally very expensive, especially if we need to sim- 
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ulate large fields of view necessary to generate mock data sets for 
next generation radio telescopes. In order to circumvent the need to 
perform such computer-intensive calculations to obtain T s , we de- 
velop an algorithm along the lines of BEARS (Tho mas et al.| 2009) 
as an approximation. In this paper we present an algorithm that 
follows the decoupling of T s from Tcmb owing to Lya photons, 
which couples the spin temperature to the colour/kinetic temper- 
ature via the Wouthuy sen-Field effect (Wouthuysen 1952). Col- 
lisional excitation and heating caused by secondary electrons re- 
sulting from hard X-ray radiation are also included. The dominant 
source of Lya flux is the background created by the redshifting of 
photons in the Lyman band into Lya. These photons are blueward 
of Lya and is injected into Lya at some distance away from the 
source. 

The amount of intrinsic Lya, ionizing and "heating" pho- 
tons is a function of the source spectral energy distribution (SED). 
Thus the evolution of the spin temperature critically depends on 
the source of reionization. Different reionization sources manifest 
themselves by influencing the IGM in markedly different ways. For 
example, deficiency of hard photons in the SEDs of "first stars", 
limit the extent to which they heat the IGM (Thomas & Zaroubi 
|2008l|Zaroubi et al.|2007l|Zaroubi & Silk|2005| ), while miniquasars 
(or miniqsos, characterized by central black hole masses less than 
a million solar), abundant in X-ray photons, cause considerable 
heating jMadau, Me iksin, & Rees|1997[|Ricotti & Ostriker| 2004 ; 
Nusser 2005 ; Furlanetto, Zaldarriaga, & Hernquist 2004 ; Wyithe 
& Loeb 2004, Furlanetto & Loeb 2002, Thomas & Zaroubi 2008, 



Zaroubi et al. 2007 ). Ionization profiles similarly have their char- 



acteristic source-dependent behavior. 

Although the question on which sources did the bulk of the 
reionization is up for debate, it is conceivable from observations of 
the local Universe up to redshifts around 6.5, that sources of reion- 
ization could have been a mixture of both stellar and quasar kinds 
(their respective roles again are uncertain). Implementing radiative 
transfer that include both ionizing and hard X-ray photons has been 
difficult and as a result most 3-D radiative transfer schemes restrict 
themselves to ionization due to stars ( [Gnedin &~A bel 2001 , |Cia^j 
|rdi et al.|[200T] |Ritzerveld, Icke, & Rijkhorst|[2003l |Susa||2006[ 
Razoumov & Cardall 2005, Nakamoto, Umemura,& Susa 2001, 
|Whalen & Norman|[2006l |Rijkhorst et al.|[2006l |Mellema et al. 
2006[ |Zahn et al.||2007[ |Mesinger & Furlanetto ||2007| |Pawlik & 



Schaye 2008). In |Ricotti & Ost riker (2004]), a 4 'semi" hybrid model 



of stars and miniqsos, like the one hinted above, has been used 
albeit in sequential order instead of a simultaneous implementa- 
tion. That is, pre-ionization due to miniqsos was invoked between 
7 ^ z ^ 20, after which, stars reionize the Universe at redshift 7. 
We in this paper would like to address the issue of simulating the 
propagation of both the UV and hard X-ray photons, exactly in 1-D 
and as approximation in 3-D. 

The focus of this paper is therefore to introduce the algorithm 
that is used to implement IGM heating in BEARS along with the 
procedure to estimate the spin temperature of the IGM. As an ap- 
plication of this technique we explore the effects of heating due to 
miniqsos, stars and, for the first time, a mixed "hybrid population". 
Subsequently, we provide quantitative and qualitative analysis of 
the differences in the 21 -cm EoR signal with and without the usual 
assumption of T s being always decoupled from Tcmb . 

The paper is organized as follows; ^describes briefly the N- 
body and 1-D radiative transfer codes used. In ^3] we describe the 
adaptation of BEARS to include Tk, followed by the calculation of 
the T s and £Tb within the simulation box. BEARS is then applied to 
three different scenarios of reionization in ^4] viz., (1) the primary 



source being stars, (2) miniqsos and (3) a hybrid population of stars 
and miniqsos. Subsequently, observational cubes of 5Tb are gen- 
erated for each of these scenarios and its properties discussed. In 
^5] we provide statistics on the simulated boxes and interpret the 
finding mainly from the point of view of the differences in 5Tb 
with and without the usual assumption that T s is decoupled from 
Tcmb- We also compare our work to that of San tos et aL| (2008 ) 
in the same section. Conclusions and discussions of the results are 
presented in ^6] along with a mention of a few topics that can be 
addressed using the data set simulated in this paper. 



2 SIMULATIONS 

In this section we briefly mention the dark-matter-only N-body 
simulations employed, the 1-D radiative transfer code developed by 
|Thomas & Zaroubi] ( |2008|> and provi de an overview of the BEARS 
algorithm designed in |Thomas et al.| ( [2"009| >. 



2.1 N-body: dark matter only 

The dark matter only N-body simulations are the same as in 
|Thomas et ah] ( [2009] ) and was performed by Joop Schaye and An- 
dreas Pawlik at the Leiden observatory. The N-body/TreePM/SPH 
code GADGET (Springel 2005 ) was used to perform a dark matter 
(DM) cosmological simulation containing 512 3 particles in a box 
of size 100 /i -1 comoving Mpc with each DM particle having a 
mass of 4.9 x 10 8 h' 1 M . 

Initial conditions for the particle's position and velocity were 
obtained from glass-like initial conditions using CMBFAST (version 
4.1; |Seljak & Za ldarriaga (1996 )) and the Zeldovich approximation 
was used to linearly evolve the particles down to redshift z = 127. 
A flat ACDM universe is assumed with the set of cosmological pa- 
rameters Q m 0.238, Q b 0.0418, fi A 0. 762, a 8 0.74, 
n s = 0.951 and h = 0.73 ( [Spergel et al.|2007| ). Snapshots were 
written out at 35 equally spaced redshifts between z — 12 and 
z = 6. Halos were identified using the Friends-of-Friends algo- 
rithm ( [Davis et al.|1985| ), with linking length 6 = 0.2. The smallest 
mass haloes we can resolve are a few times 10 10 M . The den- 
sity field was smoothed on the mesh with a Gaussian kernel with 
standard deviation oq — 2 x 100/512 comoving Mpc//i (2 grid 
cells). 



2.2 1-D radiative transfer (RT) code. 

The implementation of the 1-D radiative transfer code is modular 
and hence straightforward to include different spectra correspond- 
ing to different ionizing sources. Our 1-D code can handle X-ray 
photons, the secondary ionization and heating it causes, and as a 
result perform well both for high (quasars/miniqsos) and low en- 
ergy (stars) photons. Further details of the code can be found in 
|Thomas & Zaroubi | ( [2008] ). 

The parameter space simulated by the 1-D radiative transfer 
scheme (Thomas & Zaroubi 2008 ) include; 1) different SEDs (stars 
or miniqsos with different power-law indices), 2) duration of evolu- 
tion (depending on the life-time of the source), 3) redshifts at which 
these sources turn on and 4) clumping factor or over-density around 
a given source. Information on the ionization and temperature pro- 
files from these simulations were cataloged to be used later as in 
jThomas et al.| |2009 ). Density profiles around sources are assumed 
to be homogeneous although the code could be initialized using any 
given profile. 
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2.3 BEARS algorithm: overview 

BEARS is a fast algorithm that simulates the underlying cosmo- 
logical 21 -cm signal from the EoR. It is implemented by using the 
Nbody/SPH simulation in conjunction with a 1-D radiative trans- 
fer code; both discussed in the previous sections. The basic steps 
of the algorithm are as follows: first, a catalogue of 1-D ionization 
profiles of all atomic hydrogen and helium species and the tem- 
perature profile that surround the source is calculated for different 
types of ionizing sources with varying masses, luminosities at dif- 
ferent redshifts. Subsequently, photon rates emanating from dark 
matter haloes, identified in the N-body simulation, are calculated 
semi - analytically. Finally, given the spectrum, luminosity and the 
density around the source, a spherical ionization bubble is embed- 
ded around the source, whose radial profile is selected from the 
catalogue as generated above. In case of overlap between ioniza- 
tion bubbles, the excess ionized atoms in the overlap area is redis- 
tributed at the edge of the bubbles. For more details see Thomas 
et al. (2009). The outputs are data cubes (2-D slices along the fre- 
quency/redshift direction) of density (5), radial velocity (v r ) and 
hydrogen and helium fractions (xHi;Hii;Hei;Heii&Hein)- Each data 
cube consists of a large number of slices each representing a certain 
redshift between 6 and 11.5. Because these slices are produced to 
simulate a mock dataset for radio-interferometric experiments, they 
are uniformly spaced in frequency (therefore, not uniform in red- 
shift). Thus, the frequency resolution of the instrument dictates the 
scales over which structures in the Universe are averaged/smoothed 
along the redshift direction. The relation between frequency and 
redshift space z is given by:z = ^ — 1, where z/21 = 1420 MHz 
is the rest frequency that corresponds to the 21 -cm line. 



3 INCLUDING IGM HEATING IN BEARS 

In |Thomas et al.| p009) the BEARS algorithm was introduced; a 
special-purpose 3-D radiative transfer scheme used to simulate cos- 
mological EoR signals. In that paper we detailed the philosophy be- 
hind the code and its implementation to incorporate ionization due 
to a variety of sources. In this section, we extend the features of 
BEARS to include heating (on which the brightness temperature is 
sensitively dependent) due to hard photons within the same frame- 
work. 

The algorithm to implement heating begins in a manner sim- 
ilar to that for ionization. The Tk -profile from the 1-D radiative 
transfer code is used to embed a spherically symmetric "Tk bub- 
ble" at the locations of the centre-of-mass of dark matter haloes. 
The luminosity of the source is a function of the halo mass. The 
problem embedding a temperature profile in the simulation box is 
that, the radius of the bubble is large (> 5 Mpc) and results in ex- 
tensive overlap, even at high redshifts; a problem encountered for 
ionization bubbles only towards the end of reionization. In a quasar 
dominated part of the Universe this is particularly severe because 
of its large sphere-of-influence in heating the IGM. 

3.1 Treating overlaps: Energy conservation 

Overlap in the spheres of ionized regions was treated by redis- 
tributing the excess ionizing photons in the overlapped region uni- 
formly around the overlapping spheres ( Thomas et al. 2009). Treat- 
ing overlapped zones directly in terms of temperatures makes it dif- 
ficult to come up with any "conserved quantity". To elaborate, con- 
sider a pre-ionized zone in the simulation box close to the radiating 



source. Although at close proximity to the source, this region is 
not heated efficiently because of the absence of neutral hydrogen to 
capture the photons. On the other hand, the same source can dump 
more energy into an initially neutral zone that maybe far away from 
it. Thus, embedding the Tk -profile directly from the catalog (which 
assumes uniform IGM density) is not appropriate. 

Here we adopt a different approach where we estimate the to- 
tal energy deposited in every region of the simulation box and in- 
voke the conservation of energy deposited by the sources in the 
overlapped regions. We integrate the contribution to the energy 
budget, at a given location, from all sources whose sphere-of- 
influence extents to that location. This total energy is then redis- 
tributed equally to all the contributing sources, i.e., the energy out- 
put of each source is modified to account for the overlapping. 

To illustrate, consider Af sources that overlap at a particular 
location r and the total energy estimated at that location to be 
Etot — YliLi where Ei is the energy deposited by the i th 
source at location f. The fraction of "excess" energy attributed to 
each of the Af sources is 5E = E to t/Af. This energy 5E is then 
added to the total energy from the source to estimate a new nor- 
malization constant for the SED of the source. For example in case 
of a power-law-type source, the SED includes a normalization term 
that is determined by the total energy of the source, E to tai- How- 
ever, in the overlap case the total energy of the source is replaced 
by Etotai + SE. Now a "new" bubble of kinetic temperature is 
embedded whose profile depends on a source whose output energy 
has gotten a 5E boost. In this manner Tk(f) is estimated at every 
location in the simulation box. From this point on we need to esti- 
mate the coefficients that couple Tk to T s that will eventually lead 
to the estimation of 5Tb. It is worthwhile also to note that the re- 
sults are robust to variation in the exact algorithm used to treat the 
overlapped region. For example, if we choose to sum the energies 
deposited directly in a given region or choose max(Ei , E2 , . . . , E n ), 
where E\...E n are energies deposited by sources l...n respectively, 
the final results are very similar. 

This is obviously an approximation which assumes that the 
efficiency of the photons in heating the IGM does not depend on 
their frequency. For the case of power-law sources where most of 
the energy comes from hard photons this is probably a reasonable 
approximation. However, in the case of blackbody sources this is 
not a very good assumption because evidently the heating is done 
by a very small fraction of photons. But in these cases the overlap 
problem is not as severe as for the power-law sources since the ex- 
tent of the heating is much more limited. In summary, despite the 
approximate nature of this approach, it provides an efficient and 
reasonable resolution for the issue of the overlapping of tempera- 
ture bubbles. We note that self consistent treatment of this problem 
requires taking into account the modification of each of sources' 
SED in a different manner depending on the amplitude of its SED. 
Obviously this is computationally very expensive. 



3.2 Calculating £Tb in the volume 

Spin temperature can be thought of as a short hand notation to rep- 
resent the level population of the hyperfine states of the ground 
level of a hydrogen atom. Depending on the physical processes and 
background radiation that dominate a medium, T s is either coupled 
to the background CMB temperature or to the kinetic temperature 
of the hydrogen gas in the medium. Formally Field (1958) derived 
T s as a weighted sum of Tk and Tcmb, as; 
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1 + Vol + Vcoll 



(1) 



where y a and y co zz are parameters that reflect the coupling of T s to 
Tk via Lya excitation and collisions respectively. The efficiency of 
Lya-coupling dominates over that of collisions, especially further 
away from the source ( Thomas & Zaroubi 2008 ; Chuzhoy, Alvarez, 
|& Sh apiro 2006]). In our treatment of calculating T s , we estimate 
bothy a and y coll . 

The coefficient y co u is a function of Tk and the ionized frac- 
tion of the medium, xhii. These informations on ionization and 
heating are available from the prescription outlined in the previ- 
ous section. While simulating the 1-D profiles we simultaneously 
calculate y a as; 



Vol 
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Here, J (H) is the Lya flux density at distance r from the 
source. For miniqsos with high energy photons, Lya coupling is 
mainly caused by collisional excitation due to secondary electrons 
(Chuzhoy, Alvarez, & Shapiro 2006 ). This process is accounted for 
by the following integral, 
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where N(E;r;t) is the number of photons with energy 'E' at ra- 
dius 'r', and time 't' per unit area ( Zaroubi & Silk 2005). This num- 
ber is obtained directly from the 1-D radiative transfer by taking 
into account the the absorption due to the optical depth along the 
line-of-sight to the source at a distance 'r'. The 1-D profiles gener- 
ated are catalouged also as a function of time. Therefore the profiles 
that are embedded within the simulation are chosen to obey causal- 
ity^] The ionization cross-section of neutral hydrogen in the ground 
state is given by cr(E), fi2 = 0.416 is the oscillator strength of the 
Lya transition, e & m e are the electron's charge and mass, respec- 
tively, is the fraction of the absorbed photon energy that goes 
into excitation ( Shull & van Steenberg 1985, Dijkstra, Haiman, & 
|Loeb|2004] >. Note here that the fitting formula is a function of both 
ionization fraction and energy as given in the appendix of Dijkstra, 
|Haiman, & Loeb] ( |2004] ). 

For the case of stars, the dominant source of Lya flux results 
from the redshifting of the source spectrum blueward of Lya into 
the resonant line at different distances from the source. Frequency 
V at the redshift of emission (or the source) V is redshifted into 
VLycx atredshift z(\r\) such that; 



l + Z 

l + z(\r\)' 



(4) 



Thus, at every radius r away from the source the Lya flux 



Jo(\r\) oc N(E'-r-t)/r 2 . Where E' = 10.2 



l + z 



[eV]. 



l + z(|r|) 

Now, instead of embedding a sphere of the T s , as in the case 
of Tk, we embed instead the a bubble of Lya flux, J (|r|), as es- 
timated from Eq.[3] Since J (M) is basically the number of Lya 
photons at a given location, the overlap of two "J G bubbles" im- 
plies that the photons and hence the J Q has to be added, i.e., at a 
given spatial (pixel) location, x, y, z and time t, the total Lya flux 
is given by; 



6 Of course the catalouge has been generated with a a particular time reso- 
lution and can be made finer for more accurate results 
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(5) 



where J % (x, y, z, t) is the Lya flux contributed by i th source at the 
pixel location in the box, x, y, z and time t. Equipped with the 
quantities Tk(r), and J G (r), we can calculate y a as in Eq.[2] and 
subsequently the spin temperature through Eq. [T] Now all terms 
required for the calculation of £Tb, as in Eq.[6] are obtained. 
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4 EXAMPLE APPLICATIONS OF BEARS 

In this section we apply BEARS with its extended feature of includ- 
ing heating of the IGM, to three different scenarios of reionization. 
The models described in this section are not "template" reioniza- 
tion scenarios by any measure nor is any particular model favoured 
w.r.t the other. In fact, these models may be far from reality and 
only serve as examples of the potential of BEARS to provide a rea- 
sonable estimate of the 21 -cm brightness temperature for widely 
different scenarios of reionization. 



4.1 Heating due to stars 

The first scenario explored using BEARS is the case in which the 
sources of reionization are only stars. In this section we describe the 
model used to describe the stellar component and the prescription 
adopted to embed these stars into haloes of dark matter identified 
in an N-body simulation. 



4.1.1 Modeling stellar radiation in B EARS 

Most stars, to first order, behave as blackbodies at a given tem- 
perature, although the detailed features in the SED depends on 
more complex physical processes, age, metallicity and mass of 
the star. This blackbody nature imprints characteristic signatures 
on the IGM heating and ionization patterns ( Thomas & Zaroubi 
2008). Schaerer ('2002) showed that the temperature of the star only 
weakly depends on its mass. The blackbody temperature of the stars 
in our simulation was thus fixed at 5 x 10 4 K to perform the radia- 
tive transfer. We sample the parametre space of redshifts (12 to 6), 
density profiles around the source □and masses (10 to 1000 M©). 
The total luminosity is normalized depending on the mass of the 
star according to Table. 3 in Schaerer (2002). 

For blackbody with a temperature around 10 5 K, similar to our 
adopted value, the spectrum peaks between « 20 to ^ 24 eV 
(Wien's displacement law). Thus, from the form of the blackbody 
spectrum we can apriori expect, in the case of stars, that the ion- 
ization induced by these objects will be high, but the heating they 
cause will not be substantial given the exponential cutoff of the ra- 
diation towards higher frequencies. 



7 We can incorporate any given density profile like the homogeneous, 
isothermal or NFW profiles. For this paper however we ran all the sim- 
ulations using a homogeneous background. Overdensities around sources 
can be corrected for by using the prescribtion in |Thomas et al.| ( 2009 1. 
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Now, following the prescription in Thomas et al. (2009} we 
associate stellar spectra with dark matter halos using the following 
procedure. The global star formation rate density p*(z) as a func- 
tion of redshift was calculated using the empirical fit; 



Kinetic Temperature (log T [K]) 



P*( Z ) = P'r> 



/3exp [a(z - z m )} 



[M yr^Mpc- 3 ],^) 



1 ft — a + a exp [f3(z — a 

where a = 3/5, /3 = 14/15, z m = 5.4 marks a break redshift, 
and pm — O.15[M yr -1 Mpc -3 ] fixes the overall normalisation 
(Springel & Hernquist 2003 ). Now, if St is the time-interval be- 
tween two outputs in years, the total mass density of stars formed 



p*(z) &p*(z)5t [M© Mpc" 3 ]. 



(8) 



Notice that this approximation is valid only if the typical lifetime 
of the star is much smaller than St, which indeed is the case in our 
model because the stellar mass we input (100 M©) has a lifetime 
of about a few Myrs (Schaerer 2002), compared to a typical time 
difference of few tens of Myrs between simulation snapshots. 

Therefore, the total mass in stars within the simulation box is 
(box) ~ L\ ox p* [M©]. This mass in stars is then distributed 
among the halos weighted by their mass; 

^halo 



m*(halo) = 



M ha i Q (tot) 



M*(box) 



(9) 



where m*(halo) is the mass in stars in the dark matter "halo", 
^haio the mass of the halo and Mhaio (tot) the total mass in halos 
within the simulation box. 

We then assume that all of the mass in stars is distributed in 
stars of 100 M©, which implies that the number of stars in the halo 
is iVioo = 10 -2 x m*(halo). The luminosity of a 100 M© star is 
obtained from Fig. 1 of Schaerer (2002 ), assuming zero metallic- 
ity. The luminosity of 100 M© star thus derived is in the range of 
10 6 — 10 7 L© and this value is multiplied by iVioo t0 g et me tota l 
luminosity emanating from the "halo" and the radiative transfer is 
done by normalizing the blackbody spectrum to this value. The es- 
cape fractions of ionizing photons from early galaxies is assumed 
to be 10%. 



4.1.2 Results: Stellar sources 

Results of the evolution of the T k , T s and ST h of the IGM, when 
the sources of reionization comprise only of stars is shown in Figs 
[T] [2] and [3] respectively. 

Figure [T] shows the kinetic temperature for stellar sources at 
redshifts 10, 8, 7 and 6 at which their ionized fractions are 0.12, 
0.5, 0.83, and 0.98 respectively. The blackbody type stellar spectra 
do not have sufficient high energy photons to heat the IGM substan- 
tially far away from the source. Thus we see compact regions (bub- 
bles) of high temperature in the immediate vicinity of the source. 
In the inner parts, where the ionized fraction is high (xhii > 0.95), 
the temperatures are of the order of « 10 5 K and drops sharply in 
the transition zone to neutral IGM. The ionized region is restricted 
to about 100 kpc in physical coordinates around a halo (Tho mas &| 
|Zaroubi|2008] >. 

Figure. [2] shows T s at four different redshifts. Collisional cou- 
pling (Vcoii) is important close to the radiating source (< 200 kpc). 
The primary source of coupling between T s and Tk is the redshift- 
ing of photons blueward of the Lya line into the Lya frequency. 
It has to be emphasized here that the Lya photons are produced 
only due to the source and secondary processes. Thus we see that 
stellar sources are efficient in building up a significant background 




20 40 60 80 100 20 40 60 80 100 
Mpc/h Mpc/h 

Figure 1. Kinetic Temperature for stellar sources: Slices of Tk correspond- 
ing to redshifts 10, 8, 7 and 6 are plotted {panels left to right and top to bot- 
tom, respectively). Temperature are indicated in log T [K] on the color bar 
above. The extent of heating, both in amplitude (maximum around « 10 5 
K towards the center) and spatial extent (< 100 kpc), is extremely small. 
Basically, only the central part (ionized region) is at a high temperature and 
falls sharply during the transition into the neutral IGM. 
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Figure 2. Spin Temperature for stellar sources: Slices correspond to that 
in Fig. [T] Collisional coupling is efficient to couple the Tk to T s at the 
centre but the figures also indicate that there is sufficient J to couple Tk 
to T s well beyond the ionized region, where Tk < Tcmb- The color bar 
indicates temperature on a log T [K] scale. 



of J far away from the source which is sufficient to couple T s 
to Tk during reionization (Ciardi & Madau|20 03). As an example, 
in Fig. [12] we show the reionization history (£T D ) assuming a high 
background J Q , in other words, assuming perfect T s — Tk coupling. 

In Fig.[3](5Tb corresponding to the same redshifts as in the pre- 
vious figures is shown. In this particular model £T D values are not 
very high (\ST h \ < lOmK). At early epochs, (top-left) at z w 10, 
there are only a few sources and the heating of the IGM or the 
secondary J G is not high enough to cause substantial differential 
brightness temperatures. At later times, the ionized bubbles over- 
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Figure 3. Differential Brightness Temperature for stellar sources: Slices 
corresponds to that in Fig.[^ The temperature is now contoured on a linear 
scale with values indicated on the colorbar. In these slices, 5Tb, starting 
from zero within ionized regions attains large negative value outside be- 
cause stars are not efficient at heating the IGM beyond their ionized bubble 
and hence cannot render the IGM visible in emission. Early on, (top-left) at 
a redshift of 10, there are only a few sources and the heating induced collio- 
sional coupling or the J is not high enough to render the Universe visible, 
but the photons blueward of Lya can reach far outside the ionized region 
where they are injected into Lya frequency causing substantial coupling to 
the kinetic temperature which is much below the CMB. And, at later times, 
the ionized bubbles overlap significantly driving the brightness temperature 
to zero. 



lap significantly driving the brightness temperature to zero. Thus 
there is only a small portion of the Universe (in volume) that has 
sufficient Lya coupling and neutral hydrogen density to cause large 
differential brightness temperatures. This behaviour is specific to 
the model assumed and will be significantly different if the param- 
eters such as the star-formation rate, SED and escape fractions are 
altered. 



4.2 Heating due to miniqsos 

In this section we will consider heating due to high energy X-ray 
photons emanating from power-law type sources (eg., miniqsos). 
Because of their large mean free path (oc E 3 ) it has been diffi- 
cult to incorporate the effect of heating by power-law sources self- 
consistently in a 3-D radiative transfer simulation. 



4.2.1 Modeling miniqsos in BEARS 

Observations reveal the energy spectrum of q uasars-type sources 
typically follow a power-law of the form E~ a {Vanden Berk et al. 
2001 j JVignali, Brandt, & Schneider|2003[|Laor et al.|1997||Elvis et 
al.|1994| ). Here we assume a — 1 and thus the SED of the miniqsos 
is given by: 



1(E) = Agx E~ a 10.4 eV < E < 10 4 eV, 
with the normalization constant, 

Etotal 



Ag 



I /(tf)dE' 



(10) 



(11) 



where E to tai is the total energy output of the miniqso within the 
energy range, E range . Any complex, multi-slope spectral templates 
as in Sazonov, Ostriker, & Sunyaev (2004) can be adopted. The 
number of plausible SEDs that can be considered are numerous and 
serves as a reminder of the extent of "unexplored parameter space" 
even while considering the case of miniqsos alone, and argues for 
an extremely quick RT code like BEARS. 

The miniqsos are assumed to accrete at a constant fraction e 
(normally 10%) of the Eddington rate. Therefore, the luminosity is 
given by: 



eLedd(M) 
1.38 x 10 3 ' 



(n) 



M 



erg s 



(12) 
(13) 



The luminosity derived from the equation above is used to nor- 
malize the relation in Eq.[l0|according to Eq.[TT] The energy range 
considered is between 10.4 eV and 10 keV. Simulations were car- 
ried out for a range of masses between 10 5 and 10 9 M©. To justify 
comparing the case of stars and miniqsos we argue that, although 
the number of photons at different energies is a function of the total 
luminosity and spectral index in the case of miniqsos, if we assume 
that all photons from the miniqso are at the hydrogen ionization 
threshold, then the number of ionizing photons obtained, for the 
mass range fixed above, is about 10 50 to 10 55 ; the same order of 
magnitude as the number of ionizing photons in the case of stars 
and it also matches the numbers being employed for simulations 
by various other authors like Mellema et al. (2006) and Kuhlen & 
|Madau| ( [2005l ). 

The miniqsos are embedded into an N-body output following 
the prescription in Tho mas et al.| {2009). Dark-matter haloes are 
identified using the friends-of-friends (FoF) algorithm and masses 
of black holes assigned according to, 



M E 



10" 



M ha 



(14) 



where the factor 10 4 reflects the Magorrian relation (Magorrian 



et al. 1998) between the halo mass (Mhaio) and black hole mass 
(Mbh) multiplied with the radiative efficiency assumed to b e 10%, 
and the fraction gives the baryon ratio {Ferrarese|2002[ ). 

4.2.2 Results: miniqsos 

Contour plots of Tk, T s and and 5Th for the case of only miniq- 
sos being the sources of reionization, at four different redshifts, are 
plotted in Figs.[4][5]and[6] respectively. The high energies of the X- 
ray photons substantially affects the photo-ionization and thermal 
energy balance of the IGM. 

The marked difference between power-law type and stellar 
sources, in the extent to which they heat the IGM, is evident in 
Fig. [4] The temperatures are shown at redshifts 10, 8, 7 and 6 with 
ionized fractions 0.01, 0.54, 0.96, 0.999 respectively. In the inner 
parts of the "T k bubble" the amplitude of T k is about « 10 6 5 
K. The secondary electrons created by high energy X-ray radiation 
thermalizes to these high temperature values. Also contributing to 
extremely high temperatures in the very inner parts of the bubble 
is Compton heating ( Tho mas et al.|2009] >. The mean free path of 
X-rays are extremely large resulting in a spatial extent of heating 
that reaches many Mpcs. Thus, we see that already at redshift 10 
(top-left) the average temperature of the IGM is of the order of 10 3 
K and overlap of the temperature bubbles occur much before that of 
the ionization bubble. This example might be unrealistic, because 
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Figure 4. Kinetic Temperature for miniqso sources: Slices correspond to 
that in Fig.[T]with temperatures indicated on a log T [K] scale. Here we see 
the striking difference between miniqsos and stars in the extent of heating. 
The amplitude is about & 10 6 5 K towards the center and spatial extent of 
the heating is about a few Mpcs. Already at redshift 10 (top-left) the average 
temperature is of the order 10 3 K. 
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Figure 5. Spin Temperature for miniqso sources: Slices correspond to that 
in Fig.[T]with the temperatures indicated on a log T [K] scale. Spin temper- 
atures show an interesting ring-like behaviour (easily visible in the top-left 
panel). The reason for this is discussed in the text. 



the soft X-ray background constraint (Dijkstra, Haima ~& Loeb| 
2004) is not strictly observed, but again this serves just as an illus- 
trative example. 

Figure. [5] shows the evolution of T s as a function of redshift. 
Owing to a high secondary J Q produced due to X-ray photons, T s 
is coupled to Tk for the large part. The snapshots of T s show a 
peculiar ring-like behaviour (Refer. Fig. [5}. This is because towards 
the central part of the bubble around the source, the IGM is highly 
ionized and J Q oc xhi, obtains an extremely low value. On the other 
hand J G progressively gets lower away from the source, regions 
(>5Mpc). Therefore only a relatively narrow zone (few hundred 



Differential Brightness Temperature (mK) 




^ 60 



a 
a. 

2 40' ■ 

20 ; 

o r . . . , ....... I 

20 40 60 80 100 20 40 60 80 100 
Mpc/h Mpc/h 

Figure 6. Differential Brightness Temperature for miniqso sources: Slices 
correspond to that in Fig.[T] Note here that the brightness temperature (mK) 
in plotted on a linear scale with values as indicated on the color bar. We 
see again from the top-left panel that although the number of sources in the 
field are few, they were efficient in both increasing T^ dramatically and 
providing sufficient J Q to "light-up" the Universe in 21 -cm around them. 
The spheres in the top-left panel have 5T D of zero because they are highly 
ionized. 

kpcs) in between the ionizing front and regions further out has J 
high enough to decouple T s from Tcmb . 

Figure. [6] shows 5Tb for redshifts 10, 9, 8 and 7 in panels 
top-left to bottom-right, respectively. The 5Tb plotted in this fig- 
ure are high enough to be detectable by upcoming telescope's like 
LOFAR and MWA, whose sensitivities should reach 5Tb > 5mK 
( |Labropoulos et al.|2009) . We see in the top-left panel of this figure 
that even though the number of sources in the field is small, they 
are efficient in both increasing Tk dramatically and providing suffi- 
cient Jo to "light-up" the Universe in the 21 -cm differential bright- 
ness temperature. The 5Tb inside the spheres are zero because they 
are highly ionized (xhii > 0.99). 

4.3 Modeling co-evolution of stars & miniqsos 

In the observable Universe we know that stars and quasars do co- 
exist. Although, there are indications of the quasar number-density 
peaking around redshift two (Nusser & Silk 1993) and declining 
thereof, there are also measurements by |Fan et al.| { [2006| > of high- 
mass (> 1O 8 M ) quasars at z > 6. This allows us to envisage a 
scenario of reionization in which stars and quasars (or miniqsos) 
contributed to the ionization and heating of the IGM in a combined 
fashion. As a final example- application of BEARS, we present our 
model of a "hybrid" star-miniqso SED and examine its results on 
the ionization and heating of the IGM. 

4.3.1 The "Hybrid" SED 

The uncertainty regarding the properties and distributions of ob- 
jects during the dark ages allows us to come up with strik- 
ingly different ways of incorporating the co-existence of stars and 
miniqsos. For example, consider that there are J\f different haloes 
within our simulation box. One approach would be to suppose 
that each of these haloes host both a miniqso and stars, whose 
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Figure 7. Ionization due to the hybrid model: Slices correspond to that in 
Fig.^ Shown in these slices are the contour levels of neutral hydrogen in a 
logarithmic scale. We see here that for the hybrid model we have assumed 
the ionization proceeds quickly and the Universe is substantially ionized by 
redshift 6 (Refer Fig. [IT) for the reionization history. The extent of ioniza- 
tions interpolates between the cases of black holes and stars. 



masses/luminosity is derived according to the previous two exam- 
ples. The other approach would be to place one or two massive 
miniqsos (almost quasars), still conforming to the black hole mass 
density predictions of Volonteri, Lodato, & Natarajan (2008), and 
the rest of the haloes in the simulation box with populations of 
stars. For this example, given the minimum mass halo our simula- 
tion can resolve, we have chosen to adopt the former approach, i.e., 
a miniqso and stars in every halo. 

In our approach, given the halo mass we embed them with stel- 
lar sources as described in §4.1.1| Then the black hole mass is set as 
10 -4 times the stellar mass. Radiative transfer is performed using 
this hybrid SED, i.e., a superposition of power-law type miniqso 
SED and black body spectrum normalized according to the crite- 
rion in the previous sections. The results of ionization and heating 
are presented below. 



4.3.2 Results: Hybrid model 

The hybrid model proposed does have large amounts of high- 
energy X-ray photons from miniqsos that influence the temperature 
and ionization of the IGM prior to full reionization. The X-rays 
photons have larger mean free paths than their UV counterparts, 
and their secondary electrons have significant influence on the ther- 
mal and ionization balance. The UV photons on the other hand are 
efficient at keeping the IGM in the vicinity of the halo highly ion- 
ized. For the redshifts 10, 8, 7 and 6 we plot the ionization due to 
the hybrid model in Fig.[7]and the ionized fraction at these redshifts 
are 0.01, 0.38, 0.86 and 0.995 respectively. 

For this hybrid model the number of UV photons was not high 
enough to ionize the Universe completely by redshift 6 (the ionized 
fraction reaches 0.95). Full reionization thus could be achieved ei- 
ther by increasing the mass of the black hole and/or by changing 
the ratio between the black hole and stellar masses. 

The hybrid model considered here yielded Tk, T s and 5Tb as 
in Figs. [8] Fig.[9]and Fig.[l0] repectively. 

The hybrid model produces heating patterns that interpolates 
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Figure 8. Kinetic temperature for hybrid model: Slices correspond to that 
in Fig. [T] with temperatures indicated on a log T [K] scale. The heating 
pattern is peculiar in that the temperature shows a much "cuspier" behaviour 
towards the center. The reason is that the stellar component of the source 
ionized the very central part and the X-ray radiation can thus heat the central 
part where to even higher temperature. On the other hand the large extent 
of the heating is attributed to the power-law component of the SED. 
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Figure 9. Spin temperature for the hybrid model: Slices correspond to that 
in Fig.[T] Spin temperatures for the case of the hybrid model again shows 
the ring like behaviour as in the case of power-law sources. Temperatures 
are indicated on a log T [K] scale. 



between the results for the previous two cases and thus is qualita- 
tively different (Fig. [8}. The power-law components take over the 
role of heating an extended region and the stellar component main- 
tains a very high temperature in the central parts. 

The spin temperature (Fig. [9} again shows the characteristic 
ring-like structure around the source. The amount of secondary Lya 
radiation produced due to the miniqsos results in an efficient cou- 
pling of the spin temperature and the kinetic temperature. 

The brightness temperature for the hybrid model is plotted in 
Fig. [To] The secondary J G produced due to the power-law com- 
ponent efficiently couples the high ambient Tk of the IGM to T s 
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Figure 10. Differential Brightness temperature for the hybrid model: Slices 
correspond to that in Fig.[T] The temperatures (mK) are indicated on a linear 
scale in the color bar above. The large extent of heating and sufficient J Q 
ensures that the IGM well beyond the location of the source is rendered 
visible in 21 -cm emission. 

and this enables a large fraction of the Universe to be visible in 
the brightness temperature albeit with a few sources (top-left panel 
of the figure). The brightness temperature does not go to zero at 
redshift six because, in the model we have assumed, the Universe 
is not entirely ionized by then (for recent papers on this issue see 
|Mesinger| ( [2509t ). 

4.4 Reionization histories: A comparison 

In this section we compare the reionization histories for the three 
scenarios explored above. To create a contigious observational cube 
or "frequency cube" ( right ascension (RA) x declination (DEC) x 
redshift), we adopt the procedure described in Thomas et al. (2009). 
The procedure involves the stacking of RA and DEC slices, taken 
from individual snapshots at different redshifts (or frequency), in- 
terpolated smoothly to create a reionization history. This data cube 
is then convolved with the point spread function of the LOFAR 
telescope to simulated the mock data cube of the redshifted 21 -cm 
signal as seen by LOFAR. For further details on creating this cube 
refer to |Thomas etH1j2009) . 

As expected the signatures (both visually and in the r.m.s) of 
the three scenarios (Fig.|ll|) are markedly different. In the miniqso- 
only scenario, reionization proceeds extremely quickly and the 
Universe is almost completely (< xhii >= 0.95) reionized by 
around redshift 7. The case in which stars are the only source see 
reionization end at a redshift of 6. Also in this case, compared to 
the previous one, reionization proceeds in a rather gradual manner. 
The hybrid model, as explained previously, interpolates between 
the previous two scenarios. 

The £Tb in Fig. [TT] is calculated based on the effectivness 
of J , produced by the source, to decouple the CMB temperature 
(Tcmb) from the spin temperature (T s ). This flux, both in spa- 
tial extent and amplitude, is obviously much larger in the case of 
miniqsos compared to that of stars resulting in a markedly higher 
brightness temperatures in both the miniqso-only and hybrid mod- 
els when compared to that of the stars. However, we know that 
stars themselves produce Lya radiation in their spectrum. Apart 



from providing sufficient Lya flux to their immediate surrounding, 
this radiation builds up, as the Universe evolves, into a strong back- 
ground Jo ( Ciardi & Madau 2003), potentially filling the Universe 
with sufficient Lya photons to couple the spin temperature to the 
kinetic temperature everywhere. Thus, we plot in Fig.[l2]the same 
set of reionization histories as in Fig. [TTJbut now assuming that 
T s coupled to Tk. Clearly the "cold" regions in the Universe show 
up as regions with large negative brightness temperature. This as- 
sumption though is not strictly applicable towards the beginning 
of reionization. In ^5] we quantify the differences in the brightness 
temperature evolution when this assumption is made. 

It has to be noted that the results we are discussing here are 
extremely model dependent and any changes to the parametres can 
influence the results significantly. This on the other is the demon- 
stration of the capability and the need for a "BEARS like" algo- 
rithm to span the enormous parametre space of the astrophysical 
unknowns in reionization studies. 



4.5 Looking through LOFAR 

To test the feasibiliy of 21 -cm telescopes (specifically LOFAR) 
in distinguishing observational signatures of different sources of 
reionization, we need to propagate the cosmological 21 -cm maps 
generated in §4.4| using LOFAR's telescope response. The point 
spread function (PSF) of the LOFAR array was constructed ac- 
cording to the latest configuration of the antenna layouts for the 
LOFAR-EoR experiment. 

The LOFAR telescope will consist of up to 48 stations in The 
Netherlands|^]of which approximately 22 will be located in the core 
region (Labropoulos et al. 2009). The core marks an area of 1.7 x 
2.3 kilometres and is essentially the most important part of the tele- 
scope for an EoR experiment. The core station consists of two sets 
of antennae, the Low Band (LB A; 30-90 MHz) and the High Band 
Antenna (HBA; 1 10-240 MHz). The HBA antennae in the core sta- 
tions is further split into two "half- stations" of half the collecting 
area (35-metre diameter), separated by « 130 metres. This split 
further improves the uv-coverage. For details on the antenna layout 
and the synthesis of the antenna beam pattern refer to upcoming 
paper by Labropo ulos et al.|{2009) . 

In its current configuration, the resolution of the LOFAR-core 
is expected to be around 3 arcmins. Thus, as an example at redshift 
10, all scales below « 800 kpc will be filtered out. Fig. [T3| shows 
this effect for the reionization histories corresponding to Fig. [TT] 
The corresponding changes in the r.m.s of the brightness tempera- 
ture is also plotted. We see that, although smoothed to a large ex- 
tent, there still exists qualitative differences between the different 
scenarios. This of course is a reflection of the vastly contrasting 
models of reionization compared in this paper and difference will 
disappear if the models investigated are more similar to each an- 
other. 



5 STATISTICAL ANALYSIS OF THE SIMULATIONS 

The previous section provided a qualitative description of the dif- 
ferences in the 21 -cm £T D fluctuations for two cases i.e., with and 
without the assumption that T s = Tk. In this section, a number of 
statistical tests are performed to quantify these differences. For this 
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Figure 11. Contrasting reionization histories: From the top, reionization histories (5Tb in mK as a function of frequency or redshift) are plotted for miniqso, 
stellar and hybrid sources, respectively. The bottom panel plots the r.m.s of 5Tb as a function of redshift/frequency for all the three cases. 



purpose we focus on the scenario with the miniqsos. This exercise 
can be repeated for the other two scenarios as well. 

Firstly, Fig. [14] shows the 3-D power spectra of 5Tb at red- 
shift 12, 10, 8 and 6. In our model, reionization begins at around 
redshift 12.7 and thus z=12 can be considered the onset of reion- 
ization. The dashed curve reflects the assumption that T s = Tk 
while the solid curve takes into account the evolution of T s self- 
consistently. At the initial stages of reionization, only regions at 
close proximity to the source have been ionized and heated, and 
the volume filling by "spheres of ionization" has just begun, leav- 
ing large portions of the Universe neutral and cooler than the CMB. 
Here we assume that after the epoch of recombination the Universe 
cools adiabatically, i.e., Tk oc 1/(1 + z) 2 and Tcmb oc 1/(1 + 2?), 
which results in Tk << Tcmb at the onset of reionization. Thus, 
we have a situation in which we have assumed T s = Tk, with 
Tk << Tcmb, resulting in an artificial decoupling of T s and 
Tcmb, and is manifested as increased 21 -cm power. The discrep- 
ancy quickly (by z = 10) disappears as the "spheres of hya " 
overlap to volume-fill the simulation box and hence increase the 
T s — Tk coupling efficiency. Therefore, in our model, it is safe to 
make the assumption that T s follows Tk after about z = 10. 

From Eq. [6] it follows that 5Tb attains large negative values 
when two conditions are satisfied, T s < < Tcmb and a high neu- 
tral hydrogen overdensity. We therefore plot the mean, minimum 
and maximum values of 5Tb as a function of redshift in Fig. [15] 
As expected the mean and minimum values of 5Tb are much lower 
for the case in which the spin temperature has been set to the ki- 
netic temperature, because the Universe has cooled substantially 



below the CMB and coupling Tk everywhere in the Universe to T s 
enforces the conditions required for large negative 5Tb as stated 
above. The maximum positive values though do not differ because 
a positive 5Tb indicates that T s > Tcmb and hence has to be 
heated by the radiating source. And in regions where the source has 
heated the IGM, the presence of hya flux and colliosional coupling 
almost always guarantees the coupling of T s to Tk anyway. 

As another simple diagnostic of the efficiency of T s — 
Tk coupling we plot, (i) the differential temperature, i.e., 1 — 
Tcmb/T s , (ii) their ratio, Tk/T s , as a function of redshift (Fig. 
[T6| . Note here that if the spin temperature was artificially set to the 
kinetic temperature, the differential temperature would be saturated 
at unity. But the coupling ensures that the spin temperature remains 
below that of the CMB till about a redshift of 10. Although signifi- 
cantly different at the early phases of reionization, in the model we 
have assumed, T s catches up with Tk at a redshift of about eight, 
when the coupling due to hya and collisions become significant 
enough to tie T s to the local Tk. 

As an alternative probe to investigate the influence of hya- 
coupling coupling, in Fig [IT] we plot the power spectra, now of the 
kinetic (solid line) and the spin (dashed line) temperatures at four 
different redshifts. There are a few interesting attributes to this fig- 
ure. Firstly we see that at all scales T s is considerably lower than 
Tk at the onset of reionization (z ~ 12). This scenario changes 
rapidly as T s approaches Tk. Note that it is only towards the end 
of reionization (z=6) that T s is identical to Tk at all scales. The 
pecularity observed is that all scales in T s do not approach Tk 
simultaneously. As seen from the top two panels of the figure at 
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z — 10 and z — 8, T s approaches Tk at large scales (small k 
values). And then subsequently at small scales towards the end of 
reionization. This is because the central parts of the ionized bubble 
around the source was highly ionized and because J Q oc xhi, the 
spin temperature obtains extremely low values. Now, as you go fur- 
ther away from the source, the IGM is becomes neutral, but J G from 
the source and the secondary Lya photons is high enough to cou- 
ple the spin temperature to the kinetic temperature. Thus, the large 
scales (large volumes outside the ionized bubble) get "matched" 
first followed by the small scales, when the IGM is hot and colli- 
sions become important. 

5.1 Comparison to similar work 

Here we briefly compare the output of our simulation with that of 
Santos et al. (2008), wherein a semi-analytical scheme is developed 
to look at ionization and heating during the epoch of reionization. 
Fig[l8] shows the comparison between our simulations which has 
now been performed for the case of a power-law source starting at 
lOOeV to facilitate a balanced comparison with the lOOeV case of 
Santos et al. (2008 ). Plotted in this figure are the average gas tem- 
perature (solid line and crosses) and volume filling factor of X-ray 
radiation (dashed line and filled hexagons). The lines correspond to 
our simulation and points to Santo s et al.| ( |2008| ). Their case with 
the normalization factor f x = 1 0, P] and photon energies of 100.0 



9 | San tos et al. ( 2008) normalized their spectrum assuming a total integrated 
power density to be 3.4 x 10 40 / x ergs/s/Mpc 3 . 



eV shows similar trends. Because we start reionization a bit later 
(z — 12) the temperature is a bit lower initially. Their model is 
qualitatively different because the X-ray emissivity is tied to their 
star-formation rate. Otherwise, the manner in which heating is im- 
plemented to similar to the one described in our case albeit with- 
out terms involving compton heating/cooling. These effects though 
are important towards the central regions of the source ( [Thomas &| 
|Zaroubi|2008] ). 



6 CONCLUSIONS & OUTLOOK 

The focus of this paper was three-fold. 1) to introduce an extension 
to BEARS in order to incorporate heating (including X-ray heat- 
ing) and discuss its application to two cases of reionization and 
heating, i.e., for blackbody (stellar) and power-law (miniqso) type 
sources, 2) to monitor the evolution of spin temperature (T s ) self- 
consistently and contrast its influence against the usual assumption 
made in most reionization simulation that T s >> Tcmb, and 3) 
use the extended BEARS code to study reionization due to a hybrid 
population of stars and miniqsos as reionizing sources. 

To incorporate heating into BEARS, we embeded spheres of 
"kinetic temperature" around the sources of radiation, much like 
in the algorithm used to obtain the ionized fraction ( Thomas et al. 
2009). The BEARS algorithm implemented is extremely fast, in that 
it takes ~ 5 hours (2GHz, 16-core, 32GB shared RAM ) to per- 
form the radiative transfer (including heating) on about 35 differ- 
ent boxes (512 3 and 100 comoving Mpc) between redshifts 12 
and 6. These snapshots are interpolated between redshifts to pro- 
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Figure 13. Convolution with LOFAR - PSF: Differential brightness temperatures after convolution with a point spread function corresponding to LOFAR is 
plotted for the same case as in Fig.[n] with the temperatures plotted in mK. 



duce a contiguous data cube spanning redshifts z « 6 — 12, i.e., 
the observational range of the LOFAR-EoR experiment. The effi- 
ciency of the code facilitates the simulation of many different sce- 
narios of reionization and test their observational signatures. Apart 
from predicting the nature of the underlying cosmological signal, 
these simulations can be used in conjunction with others probes of 
reionization to enhance the detectability and/or constrain parame- 
ters concerning reionization. In Jelic et al. ( 2009} these simulations 
were cross-correlated with the simulations of the CMB anisotropy 
(mainly the kinetic and thermal SZ-effect) to probe reionization. 
In the same paper Jelic et al., repeated the same exercise by us- 
ing reionization simulations from Mellema et al. (2006 ) and found 
similar results, further emphasizing the validity of the approximate 
manner in which we solve the evolution of the kinetic temperature 
in our simulations. 

The improved BEARS code was used to simulate three dis- 
tinct reionization scenarios, i.e., miniqsos, stellar and a mixture of 
both. These simulations further emphasized that miniqsos were not 
only very efficient in increasing Tk of the IGM but the secondary 
Lya flux produced by these sources was enough to drive T s away 
from Tcmb in an extended region around the source thereby ren- 
dering the IGM around the source "bright" in 5T\>. This implies, 
if a miniqso hosting a black hole in the mass range > 10 7 M© is 
within the observing window of LOFAR, the value of £Tb and its 
spatial extent would be conducive to a possible detection. The code 
was also applied to simulate reionization with a hybrid population 
of stars and miniqsos. Every dark matter halo was embedded with 
both a quasar and a stellar component. The relation between the two 



was set according to the "Magorrian relation" we observe today. 
The total mass in quasars were constrained following black hole 
mass density estimates by Volonteri, Lodato, & Natarajan (2008). 
The effect of the hybrid population on the IGM was an interpolation 
between the scenarios of stars and miniqsos i.e., although the heat- 
ing (and ionization) was not as extended as in the case of miniqsos, 
it was definitely larger than that for stars. 

Our ultimate goal is to simulate mock data sets for the 
LOFAR-EoR experiment. Towards this end, reionization simula- 
tions performed on 35 different snapshots between redshifts 6 and 
12 were combined to form "frequency cubes" (reionization histo- 
ries, Thomas et al. (2009)) for all three scenarios i.e., stars, miniq- 
sos and the hybrid case. To fold-in the effect of the instrument we 
smoothed these frequency cubes according to the LOFAR beam- 
scale. Notwithstanding the convolution by the LOFAR beam, stark 
differences between the scenarios clearly persists. The r.m.s 5Tb 
cleary shows this difference between scenarios. 

All three scenarios of reionization discussed above was per- 
formed for two different cases, i) assuming that T s is always cou- 
pled to Tk and ii) following the evolution of T s self-consistently 
in the simulation. Figures [TT] and [12] reflect the differences be- 
tween the two cases. To further quantify the influence of the self- 
consistent evolution of T s , especially at the beginning of reioniza- 
tion, a suite of statistical analyses was performed Firstly, the 
3-D power spectra of 5Tb (Fig. [14]) for the two cases was calcu- 



10 Although these tests are specific to the scenario with miniqsos, it can 
be repeated for the case of stars and the hybrid. Miniqsos were chosen since 
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Figure 14. The 3-D power spectra of 5T h at four different redshifts (as indi- 
cated on the panels) are plotted. The solid curves indicate the power spectra 
resulting from the self-consistent inclusion of T s and the dashed curve, un- 
der the assumption T s = Tk. At the onset of reionization, z « 12 in our 
model, 5Tb is larger at all scales under this assumption because we have 
enforced, in effect, an artificial decoupling of T s from Tcmb in regions of 
the simulated Universe where Tk < < Tcmb • But by redshift 10 the power 
spectra match exactly as the Lya becomes efficient everywhere in the box 
to couple T s to Tk . Notice that the power keeps dropping as xhi — > and 
Universe gets reionized at lower redshifts. 



lated at four redshifts (12, 10, 8 and 6). We see that although there 
is a considerable difference at the onset (z=12) of reionization the 
power spectra of the brightness temperature match-up by a redshift 
of 8. This is attributed to the volume filling nature of the spin tem- 
perature coupling induced by miniqsos. The minimum, maximum 
and mean value of 5Tb as a function of redshift reveal similar dif- 
ferences between the two cases. The maximum value of 5Th is the 
same in both cases because regions where the IGM is consider- 
ably heated also corresponds to regions with a large supply of Lya 
photons and they are also regions where colliosional coupling is ef- 
fective. The minimum and average values reveal a difference since 
under the assumption that T s = Tk, regions much below Tcmb 
are artificially coupled to T s . As a final statistic we compute the 
3-D power spectra of Tk and T s and compare them (Fig. [IT] Vast 
differences do exist between the two almost till the end of reion- 
ization, when they become identical. An interesting result of this 
comparison is that the large scales of the two spectra get matched 
before the small scales in the box. The interpretation of this be- 
haviour is given in section [5] We also compared our results to that 
of |Santos et al.| ( |2008| ) and found good match for the evolution of 
gas temperature and the volume filling factor of X-ray radiation. 

In the near future calibration errors expected from a LOFAR- 
EoR observational run will be added to these simulations along 
with the effects of the ionosphere and expected radio frequency 




Redshift 

Figure 15. Plotted above are the mean (top), minimum (middle) and max- 
imum (bottom) values of 5Tb as a function of redshift. Line styles corre- 
spond to FigQ4] Because at earlier redshifts (around 12) Tk << Tcmb 
within a large fraction in the simulation box, the minimum value attained 
by assuming T s = Tk is much lower than when only regions with Lya 
excitation (which is typically near the source) is allowed to couple T s with 
Tk. And for the same reason the maximum value in both cases are similar 
because the assumption holds close to the sources as the presence of Lya 
flux ensures efficient coupling. As a result the global average 5T^ is also 
lower for the approximated model at the onset of reionization and catches 
up quickly to the "proper" model as the Universe fills up with sufficient Lya 
flux. 



interference (RFI). Subsequently, Galactic and extra-galactic fore- 
grounds modeled as in |Jelic et al.| ( [2008] ) will be merged with 
these simulations to create the final "realistic" data cube. These 
data cubes will be processed using the LOFAR signal extraction 
and calibration schemes being developed to retrieve the underly- 
ing cosmological signal in preparation for the actual LOFAR-EoR 
experiment which has already begun. 
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